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Abstract 

Computation of interlaminar stresses from the higher-order shear and normal deformable beam 
theory and the refined zigzag theory was performed using the Sine method based on Interpola- 
tion of Highest Derivative. The Sine method based on Interpolation of Highest Derivative was 
proposed as an efficient method for determining through-the-thickness variations of interlaminar 
stresses from one- and two-dimensional analysis by integration of the equilibrium equations of 
three-dimensional elasticity. However, the use of traditional equivalent single layer theories often 
results in inaccuracies near the boundaries and when the lamina have extremely large differences in 
material properties. Interlaminar stresses in symmetric cross-ply laminated beams were obtained 
by solving the higher-order shear and normal deformable beam theory and the refined zigzag the- 
ory with the Sine method based on Interpolation of Highest Derivative. Interlaminar stresses and 
bending stresses from the present approach were compared with a detailed finite element solu- 
tion obtained by ABAQUS/Standard. The results illustrate the ease with which the Sine method 
based on Interpolation of Highest Derivative can be used to obtain the through-the-thickness dis- 
tributions of interlaminar stresses from the beam theories. Moreover, the results indicate that 
the refined zigzag theory is a substantial improvement over the Timoshenko beam theory due to 
the piecewise continuous displacement field which more accurately represents interlaminar discon- 
tinuities in the strain field. The higher-order shear and normal deformable beam theory more 
accurately captures the interlaminar stresses at the ends of the beam because it allows transverse 
normal strain. However, the continuous nature of the displacement field requires a large number of 
monomial terms before the interlaminar stresses are computed as accurately as the refined zigzag 
theory. 
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1 Introdction 


Composite materials remain a topic of considerable interest to researchers due to their vulnerability 
to impact loading. Because of interlaminar bonding imperfections, delamination can be initiated by 
interlaminar stresses. This effect is amplified in sandwich composites subjected to compressive loads 
because the facesheets tend to buckle, accentuating the risk of delamination growth. The delami- 
nation can eventually result in ultimate failure of the composite. Modern composite failure criteria 
incorporate the interlaminar stresses. Therefore, accurate determination of interlaminar stresses is a 
topic of considerable interest in the research community. We will recapitulate some of the significant 
developments in this area. The state of interlaminar stress computation has been reviewed extensively 
by Kapania and Raciti [1] and Kant and Swaminathan [2]. 

Reddy [3] provided an alternative to the use of constitutive relations to obtain transverse normal 
and shear stresses by integrating the equilibrium equations of 3D elasticity directly. However, this 
technique requires in-plane derivatives of in-plane strains, a feature not explicitly available from equiv- 
alent single layer (ESL) theories such as classical lamination or first-order shear deformable theories 
in finite element implementation. For the transverse shear stress components, the equilibrium inte- 
gration approach necessitates the first derivatives of in-plane strains. For the transverse normal stress 
component, the second derivatives of in-plane strains are required. 

A significant amount of research has been conducted in the area of developing post-processing 
schemes to be used with commercial finite element software to obtain the transverse stresses. Lajczok 
[4] used a finite difference scheme to compute the higher derivatives of in-plane strains necessary 
within the equilibrium equation approach. Byun and Kapania [5] used Chebyshev and other orthogonal 
polynomials to interpolate the displacements and compute the higher derivatives of in-plane strain. Lee 
and Lee [6] developed a post-processing scheme for determining the transverse stresses in geometrically 
nonlinear composites for cross-ply, symmetric composites and sandwich composites with laminated 
facesheets and isotropic core material. Roos, Kress, and Ermanni [7] developed a post-processing 
method for determining the transverse normal stress component in doubly-curved, laminated, shell 
elements in commercial finite element codes. Noor and Malik [8] used a two-stage procedure in which 
first-order shear deformable theory (FSDT) was used to compute in-plane stresses in stage one. In 
stage two, integration of the equilibrium equations of 3D elasticity is preformed to determine through- 
the-thickness distributions of stress. An elasticity model is used to update the in-plane stresses. Stage 
two is repeated until a convergence criteria on the transverse normal stress is satisfied. Tessler and 
Riggs [9] developed a variational smoothing algorithm to accurately obtain the higher order derivatives 
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and achieve the interlaminar stresses. Tessler et al. [10] used the “Smoothing Element Analysis” for 
improving the accuracy of finite element stresses for general application and not just for interlaminar 
stresses; however, the approach is essentially a second finite element analysis. 

Recently, Slemp and Kapania [11, 12] proposed using the Sine Method based on Interpolation of 
Highest Derivative (SIHD) for problems in which the interlaminar stresses are desired. SIHD was 
developed by Li and Wu [13] as a numerical method for solving of one- and two-dimensional boundary 
value problems. The method is substantially different from traditional methods such as Ritz method 
or finite element method in which approximate the unknown function using some global or local basis 
and the basis is differentiated to obtain the derivatives. With the SIHD method, the highest derivative 
of interest is approximated using a Sine basis. The method utilizes the numerical indefinite integration 
based on the double exponential transformation proposed by Muhammad and Mori [14] to obtain the 
lower derivatives and unknown function. The concept is that by using the process of integration as 
appose to differentiation, the highest-order derivatives would be very accurate. 

Li and Wu [13] illustrated that the method is very accurate and has very good convergence proper- 
ties for one- and two-dimensional boundary value problems involving only Dirichlet boundary condi- 
tions. In [11], Slemp and Kapania assessed the accuracy of the method was examined for the analysis 
of interlaminar stresses in composite beams and plates. Numerical results were compared against 
the analytic solution for the problem of a simply-supported Timoshenko beam and for a simply- 
supported classical laminated composite plate. The method was shown to exhibit excellent accuracy 
and convergence properties for the displacements, stresses, and interlaminar stresses by using the 
through-the-thickness integration approach. 

In the next study, Slemp and Kapania used the SIHD method as a tool to solve the one-dimensional 
beam equations for the Timoshenko (FSDT) and the Bickford beam composite ESL theories. The 
numerical results for the interlaminar stresses were compared against a 3D finite element solution 
to assess the accuracy of the one-dimensional theory as compared with elasticity. The results for 
symmetric cross-ply and functionally graded composites under uniform loading indicated that the 
stresses compare well against 3D elasticity through most of the length of the beam; however, in the 
vicinity of the boundary, the results were poor. This error was attributed to the three-dimensional 
effects in the vicinity of the boundary that may not be captured by the one-dimensional ESL theories. 
The present study attempts to examine two refined beam theories in order address the viability of 
these methods as a low cost option for accurately obtaining interlaminar stresses. The present study 
also utilizes the SIHD method to solve the boundary value problem to once again illustrate the success 
with which this approach can be used to obtain the interlaminar stresses without post processing to 
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compute the necessary strain derivatives. 

There are two philosophies for improving results from ESL beam theories. First, the order of 
approximation for the displacement fields may be increased. Both the in-plane displacement and the 
transverse deflection fields could take a higher-order polynomial approximation. Such an approach is 
typically referred to as a higher-order shear and normal deformable theory (HOSNDT). HOSNDTs 
have been used by many authors to analyze composite laminates and functionally graded materials. 
The higher-order shear and normal deformable theories as derived by Batra [15] allows users the 
ability to specify the order of the assumed displacement field. Xiao et al. [16, 17] used MLPG with 
radial basis functions to analyze both thick isotropic and thick composite plates. Their results showed 
that the HOSNDT produces excellent through-the-thickness longitudinal stress variations for thick 
laminates. Transverse normal and transverse shear stresses in some cross-ply laminates computed by 
the constitutive equations were accurate using as low as a fifth-order theory. These theories have also 
been used to study functionally graded materials. Qian, Batra, and Chen [18] used HOSNDT to study 
functionally graded plates comprised of aluminum and zirconia or aluminum and a another ceramic 
such as SiC. 

A central drawback of the HOSNDT is the inability of continuous functions to approximate discon- 
tinuities. It is this philosophy that spurred the development of the so called layer- wise theories [19]. 
In layer-wise theories, the displacement is piecewise linear through the thickness. A new displacement 
degree of freedom is introduced for each layer and continuity of tractions between each layer is imposed. 
While the methods are very accurate, even comparable to three-dimensional elasticity [19, 20], the in- 
creased number of kinematic degrees-of-freedom becomes a significant drawback when a laminate has 
many layers. The “zigzag” theories attempt to introduce the piecewise linear displacement while keep- 
ing the number of degrees-of-freedom constant with increasing number of layers. Because the slope of 
the displacements are discontinuous, the strains need not be continuous through the thickness allowing 
the stresses, if desired, to satisfy the traction continuity between adhesively bonded layers. Early work 
by Di Sciuva [21] and Murakami [22] suggested displacements that satisfied the interlaminar traction 
continuity. However, the Di Sciuva theory suffers from the inability to find physical significance in the 
computed shear stress. The area integral of shear stress does not equal the shear force applied. 

Averill [23] modified Di Sciuva’s theory by introducing a penalty term in the variational principal to 
enforce interlaminar shear continuity. However, the theory suffers in much the same way as Di Sciuva. 
Averill resolves the issue at the cost of variational consistency of the boundary conditions. Recently, 
Tessler, Di Sciuva, and Gherlone [24, 20] presented a variationally consistent refined zigzag theory in 
which interlaminar shear continuity was not enforced. However, the theory serves as a substantial 
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improvement over the Timoshenko beam theory or the FSDT plate theory. Tessler, Di Sciuva, and 
Gherlone [20] compares transverse shear stresses from the constitutive relation with Pagano’s elasticity 
solution [25]. Tessler, Di Sciuva, and Gherlone also obtained interlaminar stresses by integration of 
the equations of three dimensional elasticity and showed that they very accurately approximated the 
three-dimensional elasticity solution [26]. 

In this paper, the cross-ply laminated beams are analyzed by SIHD using the HOSNDT and Tessler’s 
refined zigzag theory. The interlaminar stresses are computed using the equilibrium integration tech- 
nique and compared with a finite element solution. The remainder of this paper is arranged as follows: 
The higher-oder shear and normal deformable theory for laminated beams and the refined zigzag theory 
are presented in Sections 2 and 3 respectively with details provided for computing interlaminar stresses 
from the two theories. Two numerical examples are presented in Section 4. The paper concludes with 
a brief conclusion. 


2 Higher-Order Shear and Normal Deformable Theory for 
Beams 

The development of the higher-order shear and normal deformable beam theory (HOSNDT) following 
the general approach of Batra [15] is reviewed below. First, consider a beam with rectangular cross- 
section as shown in Fig. 1. Because the has small width ( b « a), the problem can be reduced to 
consider only the displacement field U and W. Begin by assuming a displacement field of the form, 


U{x,z) = J2z'ui(x), 
i= 0 

Tl-uj 

W(x,z ) = Y^z l Wi(x), 


(1) 


i=0 


where U and W are the longitudinal displacement and the transverse deflection respectively, and where 
the domain is x e [0, a] and z € [—h/2, h/2] . Accordingly, the strains are 


_dU_^i 

^xx o / j % 

i=0 

dW 

= -q— = )^W Z W U 


(2) 


i= 0 
n u 


dU dW i . 

lxz = -^ + -rr- = Ui + 2^z w itX . 


i = 0 


i= 0 
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Note that (i)z l 1 should be zero for all z when i = 0; however, to avoid further complication by 
defining these in a piecewise fashion, the inconsistency is noted textually. 
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Figure 1: Definition of beam geometry and coordinate system. 


The constitutive law is reduced assuming plane stress in the y direction. Note that for the traditional 
ESL beam theories studied in Slemp and Kapania [12], the stress resultants N yy , N xy , M yy , M xy , P yy , 
P xy , and Q yz were assumed to be zero where these were defined by: 


/ h/ 2 

{dyVlCxy, ZCyyi Z(X X y, Z (Tyy,Z (T X y , (Xy dZ = 0 

-h/2 


For this HOSNDT each lamina is assumed to experience plane stress (b « a) with stress components 
(jyy = a xy = tj yz = 0- Accordingly, the reduced constitutive relation is given by: 


^ XX 

°zz 


— O^e + O^e 

— ^11 t xx T t zz , 

= Ql 3 e xx + Q 33 tzz 

u xz = Q35 'Ixz 


(3) 


where is the plane stress reduced stiffness matrix for the fcth lamina (See Fig. 1). To obtain 
Qij\ in terms of the lamina engineering constants, please refer Reddy [19]. Note, the stiffnesses were 
reduced for plan stress about the y axis and not the z axis in this sense. 

The governing equations are derived using the principle of virtual work. The internal virtual work 
is written as: 

= f ^Se^dQ (4) 

J n 

where repeated index implies summation (where appropriate) and O is the beams volume. Substituting 
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Eqs. (2) and (3) into the internal virtual work and integrating by parts where appropriate yields 


SW, 


int 


b J (^11 ijUjjXx “I" ^13 ij'Mjix) bUi (TsujUj jX + ^33 ij^j) 

H” (-^51 ij'U'j “1“ -^53 ij^jix) bUi “I - -^63 ij'Mj,xx') 


dx 


”f- b {TwijUj^x -^13 ij^j) btlj ^ + b (TQUjUj TQSijWj^x') (5) 


where 


THay 


T lUj = £ qW 

k = 1 

Tllay 

T 3Hj = £ Ql3 } (*) 

fc=l 

THay 

r 6 i« = E < 5«(v) 

*:= 1 

THay 

r 8 i« = E < 5«(?) 

fc—i 


i+j+1 _ i+j+1 
Z fc+1 Z fc 


i+j + 1 


„*+j _ i+i 
Z fc+1 Z fc 




i + j 


i+j-l _ J+j - 1 
Z fc+1 Z fc 

i + j - 1 


z *+? _ z 
ffe+i z 


i+j 


i + j 


T l3ij = J2 QuU) 

fe=l 

^lat/ 

733ij = ^Q^(*j) 


fe=l 


i+j-1 _ i+j — 1 
Z fc+1 Z fc 

i+j 

i+j-1 _ i+j-1 
Z fc+1 Z fc 

i+j-1 


Tllay 


T 53ij = '£Q& ] (i) 


k = 1 

Tllay 

763ij = ]TQ$ 
fc=l 


~i+j _ ~i+j 
z fc+l z fc 

i + j 


i+j+1 i+j+1 


“fc+l 


--Zi. 


i + j + 1 


ni a y is the total number of composite layers, and repeated index of i and j implies summation over 
the terms in that index. For T 1Uj , T 13ij , T 51ij , and T 53ij , i = 0, 1, 2, ...n u ; for T 31ij , T 33ij , T 61ij , T 63ij , 
i = 0- 1 ■ 2 , . . -Ti/i ;. : for +i i ij - T 3 \ijt i I'j , and Tanj, j = 0,1,2,... //. u ; for T\ 3 iji T 33 iji 77 33 iji *md T 33 ij , 
j = 0, 1, 2, ...n TO . Also, note that if a zero-denominator is encountered in any these terms are set 
to zero. The zero-denominator only results from differentiation of a z 1 when i = 0. Therefore, these 
terms are zero and no denominator is computed. 

The external virtual work for a beam under the simply-supported condition with transverse pressure 
loading q(x) applied on the top surface of the beam (z = h/2) has the form: 

pa pa ti w 

SWext = / q(x)SW(x,h/2)dx = / —q{x) ^ [h/2) l 8wjdx 

J 0 Jo 
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The variational statement, 


b J ( T\ \ij ltj .XX T] [\ij V'j.x T T r , 1 ij Uj T ^53 ijWj^x) 8Ui 

(“t-^31 ij'Uj, T ThsijWj TQUjUj^x TQ3ijWj )X x “f" (/i/2) *) 8wi 


dx 


b ( T \ ] 'M j.x T T13 ijiXj') Siii b (TQUjUj -I - ^cY.Yij j .x ) bwi — 0, (6) 


governs equilibrium and boundary conditions. Because the variations, 8ui and Swi, are arbitrary the 
boundary value problem may be written by the following Euler equations of equilibrium: 


Siii . T] ] iij.xx T T) 3,j i^j,x T, ; : t j Vj r ['fQij'Wj,x 0 (7) 

SW{ : T31 ijUj^x T33 ijWj “h Tq 1 :jj 7 lj x “f ^(yM/j'tlYjxx = (/l/ 2 ) (8) 

For the simply support condition, the kinematics require W (0, z) = 0 and W (a, z ) = 0. This is imposed 
by constraining all Wi at the edges. At the edge x = 0, a point constraint is applied to the U degree 
of freedom ( U (0, 0) =0). This is satisfied by setting Uj(0) = 0 for i = 0, 2, 4, .... No constraints on the 
Ui for i = 1,3, 5, ... were imposed. Thus the complete set of essential and natural boundary conditions 
used for the analysis are: 


iOj(0) = 0 
Ui(0) = 0 

(7'l \ijUj x “1“ Ti3 ij'Wj') 

x—0 
Wi(a ) = 0 

(TllijUj.x T 7"i3 ijWj') 


(i — 0, 1, 2, ..., 7i w ) 

(* = 0,2,4,...) 

=0 (* = 1,3,5,...) 

(i = 0, 1 , 2 , ..., n w ') 

=0 [i = 0, 1, 2, ..., n u ) 


(9) 


where repeated index implies summation over all elements in that index. 

Note that by selecting the order of polynomial expressions for U and W, n u and n w respectively, 
the governing equations and boundary conditions are given by Eqs. (7) and (8). The boundary value 
problem may be solved using the SIHD method in an identical fashion to Slemp and Kapania [12]. 
While we omit the details, it is important to note that by using SIHD to solve the boundary value 
problems all derivatives up to Ui, xxx {x) and w i,xx are known as a by-product of the analysis and no 
additional post-processing algorithm needs to be used to find these higher order derivatives. 

The interlaminar stresses can be found by integrating the equations of 2D elasticity through the 
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thickness of the kth lamina layer (z k < z < Zk+ 1 ). 


^ - ~C{^ m ) dz+Gm 

- -III 


r ( da 


dz + 


(10) 

( 11 ) 


Note that G™, F^ k \ and H « are functions of x that may be determined from traction equilibrium 
between the layer boundaries. First, the transverse shear stresses are considered. Substituting the 
Hooke’s law and strain-displacement relations into Eq. (9), the transverse shear stress, ai z \ is expressed 
by: 

/ z I” n u n w 

Qn £ (*X**) + Qi 3 £ (iz i ~ 1 w itX ) dz + G (fc) (12) 

k L i= 0 i=l 

Performing the integration through the layer thickness, this is expressed by: 


= - 

XZ 


z i+1 - zf +1 


Qll £ ( i + l M »,xx 1 + Ql3 £ (( Z< “ 4) W i,x) 


1=0 


i = 1 


+ ai k J(x,z k ) (13) 


Differentiating with respect to x and substituting into Eq. (10) yields: 


« 2 > = - r 

J Zk 


- Qn£ 


n <* / _i+ 1 _ _i+ 1 


i—Q 


i + l 


- Ql3 £ ({z { - zi)w itXX ) + (T [k z \ x {x,Z k ) 


i = 1 


dz + H W 
(14) 


Evaluating the integration on a layer-by-layer basis through the lamina thickness yields: 


^ = Qii£ 

i=0 L 


z i+2 -4 +2 4 +1 (*-zk)' 


(i + l)(f + 2) i + l 






i= 1 


a ( x k J(x,z h )(z - z k ) + a {k J(x,z k ) (15) 


z i+1 - 4 +1 
i + l 


z l( z ~ z k) m,x 


(k)l 


Thus the transverse shear and transverse normal stresses can be found from the displacement variables. 
Theterms cr xz (a;, z k ), 0 xz , x (x, z k ), and (J zz (x, z k ) are determined by imposing continuity of the tractions 
in a layer-by-layer fashion, starting from the bottom free surface. 


3 Refined Zigzag Theory for Beam Analysis 


The development of the refined zigzag theory for beam analysis, following the approach of Tessler, 
Di Sciuva, and Gherlone [20], is reviewed below. Consider the same beam geometry shown in Fig. 1. 
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Begin by assuming a displacement field of the form, 


U{x,z) = u Q {x) + zO(x) + 4^ k \z)ip(x ), 
W(x,z) = w 0 (x), 


(16) 


where U and W are the longitudinal displacement and the transverse deflections respectively, (z) is 
a piecewise linear function of the thickness coordinate, and the domain is x 6 [0, a] and z G [—/i/2, h/2\. 
Accordingly, the strains are 


€xX = U 0 ,x(x) + z8,x(x) + 4 >( ' k \z)'ip tX (x), 
7 xz(x,z) = 8{x) + (3 (k \z)ip{x) + w QtX {x), 


(17) 


where (3^ = <j $ . Note that (3 ^ is a constant within a given layer. 

The piecewise linear variations, <f>W is set to zero on the bottom and top surfaces of the beam and 
is C° continuous at the interfaces. Therefore, the distribution over the fcth lamina may be written as: 


^ k \z) = /3^(z-z k ) + ^ k - 1) (z k ) 


(18) 


Note that (z k ) implies evaluation at z k and not multiplication. Thus determining (3^ k \ fully 

prescribes the displacement field. Rather than prescribing traction continuity between layers, Tessler, 
Di Sciuva, and Gherlone [20] requires that the following equation holds: 


qW(1 + /#)) = 


i r h/2 dz ' 

h J-h/2 Q™ 


(19) 


The assumption is made that each beam layer is in a state of plane stress in the z and y directions. 

Thus the bending stiffness, Qn\ and shear stiffnesses, , reduce to the elastic moduli, and 
(k) 

Gxz , respectively. The governing equation and boundary conditions are derived from the variational 
principle (see Tessler, Di Sciuva, and Gherlone [20]). The governing equations are: 


^xx,x A 
Afxx,x Tx “ A 

fx,X+<? = 0 


= 0 


(2A) 
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The kinematic variables are uq , 6, w, and ip. The stress resultants are defined by: 


fh/2 

{N xx ,M xx ,M <t> ,V x ,V <t> } = / (21) 

J-h/2 


To impose the simply support condition in a similar manner to the HOSNDT formulation, the kinematic 
conditions require u>o(0) = 0, wq (a) = 0, and Uq( 0) = 0 with no constraints on uo(a), ip( 0), ip (a), 9(0), 
and 9(a). Thus the complete set of essential and natural boundary conditions are: 


w 0 (0) = 0, wq (a) = 0, 
u 0 (0) = 0, N xx (a) = 0, 

M xx ( 0)=0, M xx (a) = 0, (22) 

v x (0) = 0, V x (a)=0, 

M#( 0) = 0, M+(a) = 0. 


For the present study, the SIHD method was used to approximately solve the boundary value 
problem. Again the details are omitted (see Slemp and Kapania [12]); however, it should be noted that 
doing so provides all derivatives up to uq. xxx , 9, xxx , and ip )XXX so that e xx , X x and its lower derivatives 
are known across the beam without post-processing. 

The interlaminar stresses can be derived from integration of 2D elasticity in a similar fashion to the 
HOSNDT. The following relations for transverse shear and transverse normal stresses were obtained. 


r( fc ) = -(z - Z k ) U 0 , xx ~ Ip, xx (p {k) Zk + (p {k 1} (^))] + Z 2 ^ {0’ xx + + er xz( x > z k) (23) 


,(*) 


2 

+ 


Zk(z-Zk) 


'U'Q,x 


Ip, XXX (p W Zk ~ <p {k ^(^fc))] 


z 3 - z 


k _ z k( z _ z k) ^ ^Q xxx + ip xxx p^ - a^J x (x, Zk)(z - z k ) + a^>(x,Zk) (24) 


The transverse shear and transverse normal stresses can be found from the displacement variables. 
The terms a xz (x,Zk), <J X z,x( x , z k), and a zz (x,Zk) determined from traction equilibrium between the 
layer boundaries. 
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4 Numerical solution of HOSNDT and Refined Zigzag Theory 


with SIHD 


The three-layered, [0/90/0], graphite/epoxy laminated beams of Slemp and Kapania [12] was revisited 
using the HOSNDT and the refined zigzag. This material configuration was referred to as Material A. 
The properties for this beam are give Table 1. A second beam was taken with Material B, a fictitious 
lamina with extremely different shear moduli for G 13 and G 23 . The properties are given in Table 1. 
For each case the length-to-thickness ratio was a/h = 5 where a and h are the length and thickness of 
the beam respectively. 


Table 1: Material properties for a graphite/epoxy lamina (A) and a fictitious orthotropic lamina with 
extremely different shear moduli (B). 


A 

E 1 [Msi/GPa] 24.9 / 172 

E - 2 [Msi/GPa] 1.00 / 6.89 

G 12 , G 13 [Msi / GPa] 0.50 / 3.45 
G 23 [Msi / GPa] 0.20 / 1.38 

1/12 0.25 


B 


24.9 / 172 
1.00 / 6.89 
14.5 / 100 
0.145 / 1 
0.25 


4.1 Material A 

Eight-node, linear, reduced-integration, brick elements (C3D8R) in ABAQUS/Standard were used 
to obtain an approximate three-dimensional elasticity solution for the beams. The element size was 
reduced until a similar solution from two successively fine meshes was obtained. The converged mesh 
used a total of 140 elements along the length dimension, 49 elements through the thickness and 26 
elements across the width of the beam. The elements were biased toward the ends of the beam to 
obtain a high level of mesh density near the boundaries. The simply-supported boundary condition 
was enforced at both edges by constraining all nodes on the transverse faces from U 3 translation. 
The mid-plane’s intersection with the transverse face was constrained in the Ui direction, allowing 
rotation of the simply-supported faces. The central node on a transverse edge was constrained against 
v >2 translation, preventing a rigid body mode while allowing Poisson contraction everywhere so no 
artificial loads are introduced by the boundary conditions. A uniformly distributed load was applied 
to the top surface of the beam. 

An approximate solution to the HOSNDT and the refined zigzag beam theories was obtained using 
the SIHD method as implemented in Slemp and Kapania [12]. The HOSNDT was solved using SIHD 
implemented in Fortran 90. The Sine mesh size was set by specifying the Sine span to be a = 2.0 (see 
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Slemp and Kapania [12]). A total of 81 (N = 40) Sine points were used to obtain a converged solution. 
The zigzag theory was implemented in MATLAB using 201 (N = 100) Sine points. In this case the 
Sine mesh size was also set by specifying the Sine span to be a = 2.0. 

In Fig. 2, the longitudinal and transverse normal stress components obtained by Slemp and Kapania 
[12] using the Bickford and Timoshenko beam theories, the HOSNDT with n u = 5 and 10 and n w = 3, 
and the refined zigzag theory of Tessler, Di Sciuva, and Gherlone [20] were compared with the 3D 
FEM solution. 



Bickford Timoshenko G 3D FEM Zigzag HOSNDT (n u = 5) HOSNDT (n,, = 10) 


Figure 2: Longitudinal and transverse normal stress components at x = 0.5 for simply supported 
laminated beam with Material A. The Timoshenko, Bickford, and FEM results are from Slemp and 
Kapania [12]. 

The transverse normal stress was very accurately approximated by each of the approaches and 
there is little benefit to the additional complexity of the HOSNDT or the zigzag theory. However, the 
benefit is apparent with the in-plane longitudinal stress. The Bickford theory performs very well at 
capturing the trend of the curve; however, at the material interface it is erroneous. The zigzag theory, 
while it predicts only linear through lamina distributions of in-plane stress, it matches very well with 
the FEA for the lamina interfacial stresses. The zigzag theory is superior to the Timoshenko beam 
theory and performs better than the Bickford theory even though it lacks the ability to produce the 
curved through-the-lamina-thickness stress distribution present in this thick laminate. 

For the HOSNDT, the increasing order of u displacement appears to push the interfacial longitu- 
dinal normal stress toward the FEA results; however, it is only with significant additional unknowns 
that similar accuracy as the zigzag theory is achieved. 

If Fig. 3, the transverse shear and normal stresses are plotted at the end and near the end of the 
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Figure 3: Transverse shear and normal stress components at x = 0, 0.1a for simply supported laminated 
beam with Material A. The Timoshenko, Bickford, and FEM results are from Slemp and Kapania [12]. 
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beam. In the previous study [12], it was noted that the edge effects cannot be captured by the ESL 
theories for this problem. The refined theories appear to improve substantially the transverse shear 
stress at the edge of the beam. The interfacial stress is approximated very well by the HOSNDTs and 
the zigzag theory; however, the zigzag theory is unable to capture the stress asymmetry about the 
mid-plane (higher transverse normal stress in the top lamina). The HOSNDTs appear to account for 
the stress asymmetry quite well; though neither n u = 5 nor n u = 10 (n w = 3 for both) is sufficient to 
exactly capture the stress distribution. 

While the refined theories do a better job particularly for the transverse shear stress at the edge, 
none of the present theories capture the transverse normal stress at the boundary very well. The 
HOSNDTs improve upon this stress component significantly; however, they appear to oscillate about 
the FEA solution. Further increasing the order of the in-plane and transverse deflection profiles do 
not significantly improve the results. The HOSNDT suffers from the so called Gibbs phenomenon. 
Continuous functions are never able to approximate discontinuities very well. 

At 0.1a, the refined theories each perform well at capturing both the transverse normal and trans- 
verse shear stresses. The refined theories are a substantial improvement over the ESL theories for these 
stresses. 

If the normal deformability is neglected (i.e. taking n w = 0), HOSNDT produces results very 
similar to the zigzag theory. It should be noted that the present problem is quite a difficult one in 
that the loading is uniform. The edge experiences a sudden drop in applied load. Even the 3D FEM 
solution is not plausible because the boundary conditions are not exactly satisfied. The stress should 
equal the applied load on the top surface and the derivative of the transverse normal stress with respect 
to z should be zero (by 3D equilibrium); however, neither of these are exactly satisfied by the FEA 
solution. 

4.2 Material B 

The fictitious lamina (Material B) was chosen to demonstrate a key short coming of using ESL the- 
ories to calculate interlaminar stresses. Namely, computation of interlaminar stresses from the three- 
dimensional equilibrium equations of elasticity relies on the accuracy of the in-plane derivatives of 
in-plane stresses. However, there are some lamina for which the material properties are so discon- 
tinuous that obtaining accurate derivatives of in-plane stresses is not possible. The present example 
clearly demonstrates the refined zigzag theories superiority to the Timoshenko beam theory and the 
HOSNDT. 
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Analysis was performed using SIHD with the Timoshenko beam theory, the refined zigzag theory, 
and the HOSNDT with n u = 5, 10, and 15 and with n w = 3. The Timoshenko beam theory and the 
zigzag theory were implemented in MATLAB while the HOSNDT was implemented in Fortran 90. The 
results were compared with a two-dimensional plane stress finite element analysis. For the FEA, the 
mesh size was decreased until a similar solution was obtained from two successively fine meshes. The 
converged mesh size had a total of 200 elements along the length. The elements were biased toward 
the ends with a bias ratio of 50. Through the thickness of the top and bottom layer 36 elements were 
used which were biased toward the top and bottom respectively. Though the thickness of the middle 
lamina, 12 elements were used in a uniformly distributed manner. 

The longitudinal and transverse normal stress components for material configuration B were plotted 
in Fig. 4. The bending stresses are very accurately predicted by the zigzag theory, nearly indistin- 
guishable from the FEM results. Note that the Timoshenko beam theory fails to estimate the bending 
stresses. The Timoshenko beam theory prediction for bending stress at the interface is wrong by 
about 175 % and about 40 % at the top and bottom surface. The HOSNDT with using ten monomials 
approximating the in-plane displacement performs somewhat better; however, at a substantial com- 
putational cost. With fifteen monomials, the HOSNDT results approach the accuracy of the zigzag 
theory. 



FEM - - ■ Zigzag - o - HOSNDT (n u = 5) - Q- - HOSNDT (n u = 10) - > - HOSNDT (n u = 15) 


Figure 4: Longitudinal and transverse normal stress components at x = 0.5 for simply supported 
laminated beam with Material B. 

The transverse shear and normal stresses were plotted in Fig. 5 at the edge of the beam and one- 
tenth of the length of the beam. The transverse shear stress is also very well approximated using 
the refined theories; however, the zigzag theories advantage is somewhat less convincing. While the 


16 


Thickness Coordinate, z/h Thickness Coordinate, z/h 



Figure 5: Transverse shear and normal stress components at x = 0, 0.1a for simply supported laminated 
beam with Material B. 
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interfacial shear stress is accurately approximated with the zigzag theory, the asymmetry of shear 
through the thickness cannot be captured without normal deformability. However, including normal 
deformability alone does not grant accuracy for these stress components. Note that with fifteen 
monomials approximating the U displacement, the transverse shear component at the lamina interfaces 
is still not as accurately computed using the HOSNDT than with the zigzag theory - a piecewise 
discontinuous linear theory. 

The transverse normal stress component is again inaccurate at the boundary. While introducing 
normal deformability is necessary to very accurately capture the stress at one-tenth of the length of 
the beam, the transverse normal stress at the edge is not captured very well at all by the theories 
studied. 

4.3 Summary and Conclusions 

In this paper, computation of transverse shear and transverse normal stresses or interlaminar stresses 
from refined beam theories was performed using the Sine method based on Interpolation of Highest 
Derivative (SIHD). The higher-order shear and normal deformable beam theory (HOSNDT) and the 
refined zigzag theory were introduced. Integration of the equilibrium equations of three-dimensional 
elasticity was performed to obtain the transverse normal and transverse shear stress components for 
each theory. 

The aims of this paper were to determine if the erroneous nature of interlaminar stresses computed 
by the Timoshenko and Bickford beam theories could be improved by refining the assumed displacement 
field. Like the previous study [12], the results illustrate that the SIHD method can be beneficial for 
this type of problem because the higher-order derivatives of the strain needed to perform through the 
thickness integration are accurately obtained without post processing. Two material configurations 
were considered. First, comparison was made with the material configuration of the previous study 
[12] using SIHD and the Timoshenko and Bickford beam theories. A second material configuration 
was examined which demonstrated the serious limitations of the traditional ESL theories when the 
shear moduli of lamina varies too substantially. 

The HOSNDTs were used to illustrate the effects of increasing the order of the assumed displace- 
ment field and including the transverse normal strain. However, the present implementation of the 
HOSNDT does not allow for discontinuities of the shear strains through the thickness. The numerical 
results provide very compelling evidence that the refined zigzag theory provides an excellent way to 
introduce a discontinuous displacement field without significant computational costs, and its accuracy 
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benefits substantially. The present results have shown that the refined zigzag theory provides excel- 
lent accuracy of the interfacial bending and transverse shear stresses. While the transverse normal 
stress was not improved substantially and cannot be accurately predicted without allowing a nonzero 
transverse normal strain component, the refined zigzag theory should be used over simple increasing 
the order of the in-plane displacement field. 
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